#ifndef _ERPA_H_
#define _ERPA_H_
#include"system_hf.h"
#include<eigen3/Eigen/Eigenvalues>
#include<gsl/gsl_sf_coupling.h>
using namespace std;
using namespace Eigen;


///< particle and hole twobody states, particle first
class TwoBodyState_ph: public TwoBodyState
{
 public:
  TwoBodyState_ph(int p,int h):TwoBodyState(p,h){}
  //below method is not used in rpa cal.
  bool isInAChannelWith(const TwoBodyState &) const{}
  virtual bool operator<(const TwoBodyState&) const{}
};

struct ljjtz
{
  int l,jj,tz;//jj=2*j;
  vector<int> ord_list;//a vector to store the order number of HO orbits with same l,jj and tz;
};

struct sp_info
{
  int block_number;
  int n;
};



class ERPA
{
 public:
  typedef vector<System_Tz_HF::OrbitalType> HFOrbitals;
  
  ERPA(const System_Tz_HF* Ptr);
  int set_ph_same(int J,int par,vector<TwoBodyState_ph>&ph);
  void set_AB_same(int J,int par);


  int set_ph_np(int J,int par,vector<TwoBodyState_ph>&ph);
  int set_ph_pn(int J,int par,vector<TwoBodyState_ph>&ph);
  void set_AB(int J,int par);


  bool cal_same(int J,int par);//par=0 for +, par=1 for -
  bool cal_ce(int J,int par);//par=0 for +, par=1 for -
  bool cal(int type,int J,int par);


  void set_q(int J,int par,int T);//set up onebody reduced matrix element for Q;
  double BEJ(int J,int par,int v);//B(EJT,0->v)
  double S(int J,int par,int T);
  void R(int J,int par,int T,ostream&,const double gamma=2,int points=600,const double dE=0.1);
  VectorXd q;

  
  vector<vector<TwoBodyState_ph> >ph_sames;//both the particle and hole are protons or neutrons
  vector<vector<TwoBodyState_ph> >ph_nps;//neutron particle and proton hole
  vector<vector<TwoBodyState_ph> >ph_pns;//proton particle and neutron hole

  const System_Tz_HF* pSystem_Tz_HF;
  vector<MatrixXd> A_sames;
  vector<MatrixXd> B_sames;
  vector<MatrixXd> X_sames;
  vector<MatrixXd> Y_sames;
  vector<VectorXd> Omega_sames;


  
  vector<MatrixXd> Anpnps;
  vector<MatrixXd> Apnpns;
  vector<MatrixXd> Bnppns;
  
  vector<MatrixXd> T_plus_Xs;
  vector<MatrixXd> T_plus_Ys;
  vector<VectorXd> T_plus_Omegas;

  vector<MatrixXd> T_minus_Xs;
  vector<MatrixXd> T_minus_Ys;
  vector<VectorXd> T_minus_Omegas;

  vector<MatrixXd> G_sames;
  vector<MatrixXd> G_nps;
  vector<MatrixXd> G_pns;

  int Jmax;

  void init();

  void iteration();
  double density_matrix(int bra,int ket);
  void cal_density();

  vector<ljjtz> h_ljjtzs;
  vector<ljjtz> p_ljjtzs;
  vector<sp_info> table;
  
  vector<MatrixXd> density_matrix_p;//under hf sp.
  vector<MatrixXd> density_matrix_h;//under hf sp.

  MatrixXd e_h;
  MatrixXd e_p;
  
  // private:
  double a(int bra_p,int bra_h,int ket_p,int ket_h,int J);
  double b(int bra_p,int bra_h,int ket_p,int ket_h,int J);
  double e(int bra,int ket); 
  void set_e();
};
#endif
